Flow Through Campylotic Media 
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We have found that the relation between the flow through campylotic (generically curved) media, 
consisting of randomly located curvature perturbations, and the average Ricci scalar of the system 
exhibits two distinct functional expressions (hysteresis), depending on whether the typical spatial 
extent of the curvature perturbation lies above or below the critical value maximizing the overall 
Ricci curvature. Furthermore, the flow through such systems as a function of the number of curvature 
perturbations presents a sublinear behavior for large concentrations due to the interference between 
curvature perturbations that, consequently, produces a less curved space. For the purpose of this 
study, we have developed and validated a lattice kinetic model capable of describing fluid flow in 
arbitrarily curved manifolds, which allows to deal with highly complex spaces in a very compact 
and efficient way. 

PACS numbers: 47.ll.-j, 02.40.-k, 95.30.Sf 



Many systems in Nature present either intrinsic spa- 
tial curvature, e.g. curved space, due to presence of stars 
and other interestellar media or geometric confine- 
ment constraining the degrees of freedom of particles 
moving on such media, e.g. flow on soap films [2[, solar 
photosphere Q , flow between two rotating cylinders and 
spheres to name but a few. In general, these sys- 

tems force a fluid to move along non-straight trajectories 
(curved geodesies), leading to the upsurge of non-inertial 
forces. We will denote such systems as Campylotic, from 
the greek word nafiiryXo^ for curved, media. Due to the 
arbitrary trajectories that particles through a campylotic 
medium can take, depending on the complexity of the 
curved space, the flow through these media can present 
very unusual new transport properties. Campylotic me- 
dia play a prominent role in all applications where metric 
curvature has a major impact on the flow structure and 
topology; biology, astrophysics and cosmology offering 
perhaps the most natural examples. Indeed, for several 
special cases, the flow through simple campylotic me- 
dia has already been studied, e.g. Taylor- Couette flow, 
which was originally formulated between two concentric, 
rotating cylinders [3, Q , and later extended to the case 
of spheres Q. However, beyond very simple geometries, 
the flow through more complicated structures, like ran- 
domly located stars or many biological systems, to the 
best of our knowledge, has never systematically been ad- 
dressed before on quantitative grounds. Since, in general, 
this class of flows lacks analytical solutions, their study 
is inherently dependent on the availability of appropriate 
numerical methods. Flows in complex geometries, such 
as cars or airplanes, make a time-honored mainstream of 
computational fluid dynamics (CFD), a discipline which 
has made tremendous progress for the last decades [H, Q . 




FIG. 1. Streamlines of a three-dimensional fluid moving 
through a campylotic medium. The colors denote the Ricci 
scalar R' (blue and red for low and high values, respectively). 
The gray bubbles isosurfaces stand at 1/5 of the maximum 
curvature of the system. 



However, campylotic media set a major challenge even 
to the most sophisticated CFD methods, because the 
geometrical complexity is often such to command very 
high spatial accuracy to resolve the most acute metric 
and topological features of the flow. Therefore, in this 
work, we also present a new lattice kinetic scheme that 
can handle flows in virtually arbitrary complex manifolds 
in a very natural and elegant way, by resorting to a co- 
variant formulation of the lattice Boltzmann (LB) kinetic 
equation in general coordinates. The method is validated 
quantitatively for very simple campylotic media by cal- 
culating the critical Reynolds number for the onset of 
the Taylor-Couette instability in concentric cylinders and 
spheres 0-0, El , and applied to the case of two concen- 
tric tori. 
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In this Letter, by using the new numerical scheme, we 
simulate the flow through campylotic media consisting 
of randomly distributed spatial curvature perturbations 
(see Fig. [T]). The flow is characterized by the number of 
curvature perturbations and the average Ricci scalar of 
the space. The campylotic media explored in this work 
are static, in the sense that the metric tensor and cur- 
vature are prescribed at the outset once and for all, and 
do not evolve self-consistently with the flow. The latter 
case, which is a major mainstream of current numerical 
relativity (ill 12 1, makes a very interesting subject for 



future extensions of this work. 

In order to study the campylotic media, we develop 
a lattice kinetic approach in general geometries, taking 
into account the metric tensor gij and the Christoffcl 
symbols r^.-. The former characterizes the way to mea- 
sure distances in space, while the latter is responsible for 
the non-inertial forces. The corresponding hydrodynamic 
equations can be obtained by replacing the partial deriva- 
tives by covariant ones, in both, the mass continuity and 
the momentum conservation equations. After some alge- 
braic manipulations, the hydrodynamic equations read as 
follows: dtp + (pu l )-i — 0, and dt(pu l ) + T l j — 0, where 
the notation -i denotes the covariant derivative with re- 
spect to spatial component i (further details are given in 
the Supplementary Material [X3j ) . The energy tensor T u 
is given by, T u = Pg %3 + pu % v? — p{g l3 u\-\- g u 3 .,+ g lJ u\), 
where P is the hydrostatic pressure, u l the i-th con- 
travariant component of the velocity, g l ° the inverse of 
the metric tensor, p is the density of the fluid, and fi is 
the dynamic shear viscosity. 

Since lattice Boltzmann methods are based on kinetic 
theory, we construct our model by writing the Maxwell- 
Boltzmann distribution and the Boltzmann equation in 
general geometries. The former takes the form (lij ]: 
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where g is the determinant of the metric g^ , and 9 is the 
normalized temperature. The macroscopic and micro- 
scopic velocities, u l and are both normalized with the 
speed of sound c s = -^//cbTo/to, feg being the Boltzmann 
constant, To the typical temperature, and m the mass of 
the particles. Note that the metric tensor appears ex- 
plicitly in the distribution function, due to the fact that 
the kinetic energy is a quadratic function of the velocity, 
u l Ui = g,ijU l ui . To recover the macroscopic fluid dynamic 
equations, we have to extract the moments from the equi- 
librium distribution function. The four first moments of 
the Maxwellian distribution function on a manifold are 
given by, 



p = / f<% , P u* = / f?dz 



pu U J 



pOg 1 
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These moments are sufficient to reproduce the mass and 
the momentum conservation equations. Here, for simplic- 
ity we have used <i£ to denote d^d^dt; 3 and the Jaco- 
bian of the integration is already included in the Maxwell 
Boltzmann distribution, through the determinant term 

In the absence of external forces, in the standard the- 
ory of the Boltzmann equation, the single particle distri- 
bution function /(x l , £ l , t) evolves, according to the equa- 
tion, dtf + C^if = C(f), where C is the collision term, 
which, using the BGK approximation, can be written as, 
C = — (1/t)(/ — / cq ), with the single relaxation time r. 
This equation can be obtained from a more general ex- 
pression, df/dt = C(/), where the total time derivative 
now includes a streaming term in velocity space due to 
external forces, f = d t f + ^-dj + %-d p if, with p l the 
i-th contravariant component of the momentum of the 
particles. Using the definition of velocity, = dx t /dt, 
and due to the fact that the particles in our fluid move 
along geodesies, which implies the equation of motion 



we can write the Boltzmann equation as [l5| 



(3) 



(4) 



where we have used the definition of the momentum, 
p l = rn^ 1 . Note that the third term of the left hand side 
carries all the information on non-inertial forces. Thus, 
all the ingredients required to model a fluid in general ge- 
ometries within the Boltzmann equation are now in place. 
Note that the Christoffcl symbols and metric tensor are 
arbitrary and therefore we can model the fluid flow in 
curved spaces, whose metric tensor is very complicated 
and/or only known numerically. 

Since the contravariant components of the velocity are 
free of space-dependent metric factors, they lend them- 
selves to standard lattice Boltzmann discretization of ve- 
locity space. All the metric and non-inertial information 
is conveyed into the generalized local equilibria and forc- 
ing term, respectively. These features are key to the LB 
formulation in general manifolds. As an additional fea- 
ture, complex boundary conditions related to a specific 
geometry, e.g. surface of sphere, in many cases, can be 
treated exactly by cubic cells in the contravariant coor- 
dinate frame, thereby avoiding stair-case approximations 
typical of cartesian grids. The details of the discretiza- 
tion of this model on a lattice can be found in the Sup- 
plementary Material 13 [. 

To provide numerical validation of our model, we study 
the flow through one of the simplest campylotic medium, 
the Taylor-Coutte instability in three different geome- 
tries, i.e. two concentric rotating cylinders, spheres and 
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FIG. 2. Critical Reynolds number Re c , as a function of the 
parameter r\ — a/b at the onset of the Taylor-Couette in- 
stability, for two concentric rotating cylinders (red) and tori 
(blue). Theoretical values for the case of the cylinders agree 
with Ref. The left inset shows the critical Reynolds num- 
ber for the case of two concentric spheres, and the two colored 
spheres the radial and axial components of the fluid velocity 
for the spherical case. Blue and red colors denote low and 
high values, respectively. 




FIG. 3. Flux reduction $o — with respect to the flat case, 
as a function of the number of curvature perturbations for 
ao — 0.01 and ro = 2.0. The solid line is the analytical curve 
according to Eq. ([5]). Shown in the inset is the normalized 
average curvature scalar of the space, R/Rmax, and the nor- 
malized reduced flux 1 — <I>o as a function of ro . Both Ricci 
scalar and flux reduction exhibit a maximum at intermediate 
values of ro. Since the two maxima are slightly shifted with 
respect to each other, the reduced flow as a function of R 
exhibits an hysteresis loop (see next Figure 2]). 



tori, respectively. Full details of the validation are given 
in the Supplementary Material[13j. In Figure [21 we re- 
port the critical Reynolds number as a function of the 
aspect ratio rj = a/b, where a and b are the minor and 
major radii, respectively. As one can appreciate, for the 
cylindrical geometry we obtain excellent agreement with 
analytical theory [5(, and a similar match with experi- 
mental data Q is found for the spherical case. We have 
also computed the torque coefficient, and found reason- 
able agreement, within a few percent, with experimental 
data [lfl, [l7f. For the case of two concentric rotating 



tori, the critical Reynolds numbers for different config- 
urations can also be observed in Fig. [21 showing values 
around 10% larger than for the case of cylinders. Further 
details can be found in the Supplementary Material [l3| . 

Next, we move to a genuinely campylotic medium, con- 
sisting of randomly located curvature perturbations. To 
this purpose, we define a coordinate system (x, y, z), such 
that its metric tensor takes the form: gij = #y(l — 

a o 2n=i cx P(~ r «/ r o))i where n labels each local curva- 
ture perturbation located at r n = (x n ,y n ,z n ), N is the 
total number of perturbations, r n = \r n \, and ro char- 
acterizes the size of the deformation. Note that the co- 
efficient ao can be either signed, depending on whether 
a positive or negative curvature is imposed, respectively. 
In our study, we have chosen to work with positive val- 
ues of ao, due to the analogy with a system of randomly 
located stars, which produce deformations in the met- 
ric tensor of spacetime[l[. The Christoffel symbols are 




FIG. 4. Flux reduction, $o — as a function of the aver- 
age curvature, R, for large and small values of ro. We have 
fixed a = 0.00002 and N = 1024. The solid lines denote the 
analytical curves according to Eqs. ^ and J7J. The inset 
shows the hysteresis loop which arises by parametrizing the 
flux-curvature relation in terms of ro. Here, r c is the radius at 
which the Ricci curvature attains its maximum upon increas- 
ing ro. The lower and upper branches correspond to r c < ro 
and r c > ro, respectively. 
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calculated numerically. The flux is calculated by the ge- 
ometrical relation, $ = J s pu x ^g xx gdS, where S is the 
cross section at the location where the measurements are 
taken. Since the fluid dynamic equations only contain the 
metric tensor and its first derivatives (via the Christoffcl 
symbol), and due to the fact that particles move along 
geodesies according to Eq. ([3]) , it is natural to expect that 
the flow could be characterized by a quantity that con- 
tains the metric tensor and its first derivatives. Although 
the Christoffel symbols T % - k meet this requirement, they 
are not components of a tensor, and therefore they are 
not invariant under a coordinate system transformation 
(physics should not depend on the choice of the coor- 
dinate system). An invariant, or tensor, that can be 
used to characterize the system is the Ricci tensor Rn-. 
In this work, we use the Ricci scalar (curvature scalar) 
which can be calculated from the Ricci tensor, Rij, by 
contraction of the indices, R' = g iJ Rij. The relation be- 
tween the metric tensor and Christoffel symbols and the 
Ricci tensor can be found in the Supplementary Material 
fl3j . To study this particular system, we use a lattice 
size L x x L y x L z of 128 x 64 x 64, and r = 1. All 
quantities will be expressed in numerical units. To drive 
the fluid through the medium, we add an external force 
along the x-componcnt, which in all simulations takes the 
value, f ex t = 5 x 10~ 5 . The flux in flat space, i.e. in the 
absence of curvature perturbations is denoted by <E>o- 

Shown in Fig. [H are the velocity streamlines, the Ricci 
scalar R' and the high-curvature locations, represented 
by gray isosurfaces. Note that the streamlines are very 
complex, as the flow can orbit around the spheres before 
continuing its trajectory [3, [3]. Also we can see how the 
curvature perturbations interact, creating non-spherical 
shaped isosurfaces. 

Fig. [3] shows the flux reduction $o — ®, as function of 
the number of curvature perturbations, N. We observe 
that the flux $ decreases with N. This effect is due to the 
interplay between the longer trajectories that particles 
must take and the acceleration due to the non-inertial 
forces, see Eq. (J3j) . Note that, in general, for systems 
with different configurations (e.g. negative ao), we could 
expect that the combination of the two effects might lead 
to higher flux by increasing N. We also see that the flux 
depends linearly on N for low concentration of curvature 
perturbations, and only sublinearly at higher concentra- 
tions. This is due to the fact that at low concentration, 
the average distance between curvature perturbations is 
large, and consequently each perturbation adds up as a 
single modification to the total spatial curvature. How- 
ever, as the concentration is increased, the curvature per- 
turbations start to interfere with each other and conse- 
quently the space becomes less curved (decrease of the 
overall Ricci curvature). The flux is found to obey the 



following law, 



$ - $ = A\ 



N/N 



1 + (N/N ) 2 



(5) 



where A\ = 163 ± 2 and 7V = (1.54 ± 0.03) x 10 4 are 
fitting parameters. The parameter Nq denotes a charac- 
teristic number of curvature perturbations, above which 
the sublinear behavior sets in (N > No). 

In the inset of Fig. [3J we observe that by fixing the 
number of curvature perturbations N = 1024 and the 
strength ao = 0.01, and changing the range of the per- 
turbation, ro, the difference $o — $ presents a maxi- 
mum for a given ro ~ r c . Furthermore, another inter- 
esting result is that the average curvature, here defined 
as R = — 10 8 < R' > (where < ... > means average 
over space), shows the same qualitative behavior. Since 
by increasing ro the metric tensor components decrease 
monotonically, this maximum is due to the Christoffel 
symbols (or non-inertial forces), which can be character- 
ized via R. However, the maxima are slightly shifted, 
due to the fact that the Ricci scalar does not uniquely 
determine the metric tensor and Christoffel symbols, the 
quantities that play a key role in the fluid dynamic equa- 
tions. Taking into account this effect, we can plot the 
flux reduction $o ~ ^ as a function of R, and find that, 
indeed, for ro < r c , the flux decreases by increasing the 
average curvature R with a different law than for the case 
of large values of ro > r c (see inset of Fig. [3} . This gives 
rise to a hysteresis-shaped curve, the reason for this hys- 
teresis being that the metric tensor is different for ro < r c 
and ro > r c , even if R takes the same value. However, 
in both cases, the system shows that higher values of the 
average curvature R always result in a lower flux. The 
behavior of the flux for ro < r c is well represented by the 
following law: 



* * A R (, R 

$o-* = A 2 — 1 + — 

ito V -^0 



and for the case of ro > r c 



$o - <1> = A : - 



R 



Ro + R 



(6) 



(7) 



where R = 5.2 ± 0.1, A 2 = 50 ± 2, A 3 = 154 ± 4, and 
<£' = 5 ± 1. The quantity Rq is related to the maximum 
curvature achieved by the system and the intersection of 
the two laws (see Fig. [4j) . The other interesting quantity 
is which represents the difference of flux between ro 
r c and ro <C r c , when the curvature scalar becomes zero, 
and it is due to the fact that in both cases, although 
the space has no curvature, it has nonetheless different 
metric tensors. 

Summarizing, we have explored the laws that rule the 
flow through campylotic media consisting of randomly 
distributed curvature perturbations, and shown that, for 
the configurations studied in this Letter, curved spaces 
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invariably support less flux than flat spaces. Further- 
more, the flux can be characterized by the Ricci scalar, 
a geometrical invariant that contains the metric tensor 
and Christoffel symbols, the quantities that appear in 
the fluid dynamics equations. The trajectories of the flow 
can become very complicated due to the total curvature 
of the medium, presenting, in some cases, orbits winding 
several times around regions with high curvature. The 
present method opens the possibility to apply the actual 
model to astrophysical systems, where the curvature of 
space is due to the presence of stars and other interstel- 
lar material. We have not considered time curvature, 
since its contribution remains sub-dominant unless mass 
is made extremely large. 

To calculate the flux in campylotic media, we have de- 
veloped a new lattice Boltzmann model to simulate fluid 
dynamics in general non-cartesian manifolds. The model 
has been successfully validated on the Taylor-Couette 
instability for the case of two concentric cylinders and 
spheres, the inner rotating with a given speed and the 
outer being fixed. We also studied the Taylor-Couette in- 
stability in two concentric rotating tori, finding that the 
critical Reynolds number for the onset of the instability 
is about ten percent larger than the one for the cylinder. 
By solving the Navier-Stokes equations in contravariant 
coordinates, which can be represented on a cubic lattice 
precisely in the format requested by the lattice Boltz- 
mann formulation, the present model opens up the pos- 
sibility to study fluid dynamics in complex manifolds by 
retaining the outstanding simplicity and computational 
efficiency of the standard lattice Boltzmann method in 
cartesian coordinates. The case of dynamically adaptive 
campylotic media, in which the metric tensor and curva- 
ture would evolve self-consistently together with the flow, 
makes a very interesting subject for future extensions of 
the present lattice kinetic method in the direction of nu- 
merical relativity 2(], 2l| . 

The authors are grateful for the financial support of the 
Eidgcndssischc Technischc Hochschulc Zurich (ETHZ) 
under Grant No. 06 11-1. 



Covariant derivative and Ricci Tensor 

The formulation of fluid equations in general coordi- 
nates implies the replacement of partial derivatives with 
the corresponding covariant ones. Given a vector A 1 , the 
covariant derivative is defined by 



d 3 X 



Y) k A k 



(8) 



where r* fe is the Christoffel symbol associated with 
the curvature of the metric manifold, namely Y l j k = 
y m {^Br + ^ - For an arbitrary tensor of 



second order, the covariant derivative is given by 

Af = diA tk + Y l ml A mk + Y k ml A im . (9) 



Here and throughout, according to Einstein's convention, 
repeated indices are summed upon. 

The Ricci tensor Rik is related with the metric tensor 
and Christoffel symbols by the relation, 



r)Y l 



dx l dx h 



1 ik L Ir 



l 1 kr. 



(10) 



Tensor Hermite Polynomials 

The Lattice Boltzmann formulation in general geome- 
tries makes strong reliance on Hermite expansion of the 
kinetic distribution function. The first three Hermite 
polynomials arc, 



H {Q) - 1 



(lla) 



(lib) 



flic) 



H l ^ k = crC - (S 1J £, k + S kj C + 5 lk ^) , (lid) 
where we have used the Kronccker delta 5 lJ . 



Supplementary Material 

We show the details of the new lattice kinetic model to 
study campylotic media, and include a respective valida- 
tion by studying the Taylor-Couette instability for the 
case of two concentric rotating cylinders, spheres and 
tori. We also implement a convergence study showing 
that the model presents nearly second order convergence, 
and introduce basic relations in differential geometry like 
the calculation of covariant derivatives and the Ricci ten- 
sor. 



Hermite polynomials expansion 

Let us expand the distribution function f (x l , ^ , t) in 
the form, 

f(x i ,e,t)=w(oJ2^ a (n)(^,t) H (n ) (e) , (12) 

n=0 

where the coefficients ai n ) are n-th order space-time de- 
pendent tensors, and Ht n ) arc the tensorial Hermite poly- 
nomials of n-th order. The weights w(£) are defined as: 

^o-j^r^M-em ■ (13) 
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The coefficients a(„) can be calculated with the relation, 



(14) 



(„) = / fH( n )(ZH 



To recover the correct hydrodynamic equations, the 
model must be built in such a way as to recover the first 
four moments, given by, 



P 



p9(u i g* k + u j g ik + u k g ij ) + pu i u j u k 



(15a) 



(15b) 



(15c) 

The fourth one ensures that the dissipation term 
achieves the correct form. To this purpose, we need to 
expand the distribution function at least up to the third 



order Hermite polynomial (The explicit expression of the 
Hermite polynomials have been given above). Thus, us- 
ing Eq. (|14[l. and replacing the Maxwell-Boltzmann dis- 
tribution for a manifold, we obtain: 



f cq_ V3P 
J ,~ -3/2 eXP 



(2nd) 



1 

26 



. (16) 



Next, by taking 9 = 1 (isothermal limit), we obtain the 
coefficients of the expansion, as follows: 



f(o) = P J a (i) = P u ! 



',,J 



pu u 



(17a) 



(17b) 



aja) = (s ?J - 513 > fe + (g kj S kj )u l + {g lk ~6 lk )u j +pu i u j u k . 

(17c) 

Therefore, the truncated equilibrium distribution func- 
tion up to third order, using Eq. (TT2"j) . reads as follows: 



J 



.(§ + 2 ^v + \?9 ij ? - \ee + \(?rf) 2 - \? 



2 6^ ; 



(18) 



With this, we have expanded the equilibrium distribution 
function up to third order in Hermite polynomials. Next, 
we need to expand also the forcing term, Vj^^d^f, in 
the Boltzmann equation, 

dtf + €d l I-T) k ^ed e f = C(f) , (19) 

Due to the fact that the distribution function can be 
written using Eq. (|12|) . and invoking the properties of 

I 



I 

the Hermite polynomials, 

wH\ n) = (-\) n d?w , (20) 
we can write the forcing term as, 

00 a F i 

F %f = W H^W H l) . (21) 
n=l v 

where we have introduced the notation, F % = — r*. fe £ J £ . 
Then, replacing the coefficients from Eq. (fl7|) . and the 
corresponding Hermite polynomials, Eq. (|11[) . we obtain 
the forcing term, 



(22) 

I 

With every term expressed as a series of Hermite poly- tion according to standard Hcrmite-Gauss projection of 
normals, all is in place to proceed with the LB discrctiza- the continuum Boltzmann equation. 
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Lattice Discretization 

In order to formulate a corresponding lattice Boltz- 
mann model, we implement an expansion of the Maxwell- 
Boltzmann distribution in Hermite polynomials, so as to 
recover the moments of the distribution function up to 
third order in velocities, as it is needed to correctly re- 
produce the dissipation term in the hydrodynamic equa- 
tions. The expansion of the Maxwell-Boltzmann distri- 
bution was introduced by Grad in his 13 moment system 
[iil ] . Since this expansion is performed in velocity space, 
and the metric only depends on the spatial coordinates, 
we expect such an expansion to preserve its validity also 
in the case of a general manifold. We have followed a 
similar procedure as the one described in Refs. 23, III- 

For the discretization of the Maxwell Boltzmann distri- 
bution ([To]) and the Boltzmann equation (fTOj) . we need a 
discrete velocity configuration supporting the expansion 
up to third order in Hermite polynomials. Our scheme 
is based on the D3Q41 lattice proposed in Ref. fijj ]. 
which corresponds to the minimum configuration sup- 
porting third-order isotropy in three spatial dimensions, 
along with a H-theorem for future entropic extensions 
[26j of the present work. 

In the following, we shall use the notation c\ to de- 
note the z-th contravariant component of the vector num- 
bered A. Thus, the discrete Boltzmann equation for our 
model takes the form, f\{x l + c\St,t + St) — f\(x l ,t) = 
— ^fifx — f^+btFx, where T\ is the forcing term, which 
contains the Christoffcl symbols, and is the discrete 
form of the Maxwell-Boltzmann distribution, Eq. (p~6]l . 
The relevant physical information about the fluid and 
the geometry of the system is contained in these two 
terms. The macroscopic variables are obtained accord- 
ing to the relations, p = S a L /a, pu l = J^xLof^l- 
The shear viscosity of the fluid can also be calculated as 
M = P( T ~ l/2)c2<Jt. 

In the following, we shall use the notation c\ to denote 
the vector number A and the contravariant component 
i. The cell configuration D3Q41 has the discrete veloc- 
ity vectors: (0,0,0), (±1,0,0), (±1,±1,0), (±1,±1,±1), 
(±3,0,0), (0,±3,0), (0,0, ±3), and (±3, ±3, ±3). The 
speed of sound for this configuration is c 2 s = 1 — \/2/5. 
With this setup, and taking into account that the vec- 
tors £* and u l are normalized by the speed of sound, we 
obtain the following equilibrium distribution, 



jreq 
J A 



2 ci 



1 ( c >y _ 1 u _ w 1 (c>y 

2 c A s 2 9 2 c 2 s 6 c 6 s 

1 (c>')(^V) 1 (cjy) k k 

9 ~i h 9 „4 ^ c a5 c a 



(23) 



2 c; 
1 (c>*) 



(g jj - 3) 



u-g 



where the weights w\ are defined as, W(o.o.o) = 
^(5045 - 1507^10), W(1)0 ,o) = ^5 - i, u>aii l0) = 
i(55-17^T0), ^(1,1,1) = ^(233^10-730), w (3 ,o,o) = 
^(295 - 92^10), and 10(3,3,3) = I5 £ 55 (130 - 4lVlO). 

The discrete Boltzmann equation for our model takes 
the form, 



fx (x* + c\St , t + St ) - f x (x { , t ) = - - (/ A - f c x q ) 

T 

where we have introduced the forcing term, 

(c[u l )c\F{ u>F{ 



StF x , 
(24) 



StF\ = w\p 





- S kl 4 


k I 

u u 




An 


c™ 




4 


c 



c k xc[c\F{ 



(25) 



with F{ = — rj fc £^£* and S kl is the Kronecker delta. In 
the presence of an external force F ext , this simply extends 
to Fl^Fl + F^ tx . 

In order to recover the correct macroscopic fluid equa- 
tions, via a Chapman-Enskog expansion, the other mo- 
ments, Eq. (TT5]), also need to be reproduced. A straight- 
forward calculation shows that the equilibrium distri- 
bution function f x q meets the requirement. The shear 
viscosity of the fluid can also be calculated as p = 
p(r — l/2)c 2 s 8t. In this way one can calculate the fluid 
motion in spaces with arbitrary local curvatures. 



Convergence Study 

To check the convergence of the model, we simulate the 
Poiseuillc profile for the velocity on a two-dimensional 
ring. For this purpose, we use the metric tensor in polar 
coordinates, g rr = 1, ggg = r 2 , and g zz = 1, where r 
is the radial coordinate, 9 is the azimuthal angle, and z 
the axial coordinate. Thus, the non- vanishing Christoffcl 
symbols for this metric are given by, T r ge = —r, and T e rS = 

a = l/r. 

Our system consists of a two-dimensional ring with in- 
ner radius a and outer one b. On this ring, we impose a 
constant force f a in the ^-direction. For the simulation 
we choose r = 0.6. The forcing term f a is set to 0.05. 
All numbers are expressed in numerical units. The in- 
ner radius of the ring is taken as a = 1.0 and the outer 
radius as b = 1.064. We have taken periodic boundary 
conditions in the direction 9 and z, and free boundary 
conditions at r = 1.0 and r = 1.064. 

To obtain a quantitative measure of the conv erg ence 
we use the Richardson extrapolation method [27], ■ In 
this method, given any quantity A(Sx) that depends on 
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FIG. 5. Relative convergence error as a function of the num- 
ber of grid points. Here, the relative error is calculated by 
taken the mean value of the relative errors at every location 
grid point. 



a size step Sx, we can make an estimation of order n of 
the exact solution A by using 



on A ( Sx 

A = lim A(Sx) « 

5x-K) v ' 2" - 1 



A(6x) 



+ 0(5x n+L ) 



(26) 

with errors 0(5x n+1 ) of order n + 1. Thus the relative 
error between the value A(Sx) and the "exact" solution 
A can be calculated by 



E r (5x) 



A(Sx) - A 



A 



(27) 



In our case, the quantity A is the fluid density p, when 
the fluid reaches the steady state, and we set up n = 2. 
Indeed, the relative error with respect to the "exact so- 
lution" decreases rapidly with increasing grid resolution 
(see Fig. [SJ and we can see that the present scheme ex- 
hibits a near second-order convergence. This is basically 
in line with the convergence properties of classical LB 
schemes. 



Validation 

To provide numerical validation of our model we study 
the Taylor-Couette instability, which develops between 
two concentric rotating cylinders. We calculate the criti- 
cal Reynolds number, Re c , which characterizes the tran- 
sition between stable Couette flow and Taylor vortex 
flow. To this purpose, we use the metric tensor for cylin- 
drical coordinates (r, 9, z), g rr = 1, ggg = r 2 , and g zz = 1, 
where r is the radial coordinate, 9 is the azimuthal an- 
gle, and z the axial coordinate. Thus, the non-vanishing 



Christoffcl symbols for this metric are given by T 7 eg = —r, 
and T% = Tl = 1/r. 

In our system, the inner cylinder has radius a and the 
outer one radius b. We performed several simulations, by 
varying the Reynolds number for different aspect ratios 
rj = a/b. The Reynolds number, assuming that the outer 
cylinder is fixed, can be defined as Re = (aS / v)d6 / dt 
where dO/dt is the angular speed of the inner cylinder 
and 5 = b — a. The inner radius a is always set to a = 1, 
and for a given value of rj, the outer radius b and 5 are 
calculated. In order to vary Re, at fixed rj, we change 
the angular velocity of the inner cylinder. For this sim- 
ulation, we use a rectangular lattice of 128 x 1 x 256 
cells and choose r = 1 (all values are given in numerical 
units). We use periodic boundary conditions in the 9 and 
z coordinates. At r = a and r = b boundaries, we have 
used free boundary conditions, together with a condition 
to impose the respective angular velocity at each bound- 
ary by evaluating the equilibrium function with those 
values. Note that the boundary conditions can be imple- 
mented as if they referred to a cartesian geometry, due 
to the use of contravariant coordinates, leading to an 
approximation-free representation of curved geometries. 
For smooth manifolds, the new scheme is about three 
times slower than a standard cartesian version, which is 
mainly due to the calculation of the metric and curvature 
terms, as well as to the use of third order equilibria to 
enhance stability. Clearly, the advantage of the present 
scheme lies in the treatment of complex manifolds which 
would require very high cartesian grid resolution. 

In Fig. [2 we can observe the critical Reynolds num- 
ber as a function of rj, as predicted by the simulation 
and compared with the theoretical values from Ref. [H, 
finding excellent agreement. We have implemented the 
same simulation using a lattice size 64 x 1 x 256 cells in 
order to study the influence of the boundary conditions, 
and we found an error of around 2.5% respect to the 
theoretical values, which is a clear evidence of the sen- 
sitivity to the boundary condition implementation. We 
have been able to simulate Reynolds number of around 
7000 by using r = 0.55. Note that our model works in 
contravariant coordinates and due to the presence of a 
metric tensor, the time step is not necessarily unity. For 
this reason, even when the relaxation time is not small, 
the computed kinematic viscosity can achieve very small 
values, leading to large Reynolds numbers. 

For the case of two rotating spheres, we consider the 
inner sphere with radius a and the outer one with ra- 
dius b. We use standard spherical coordinates (r,<j>,6), 
being r the radial, 4> the azimuthal, and 9 the polar co- 
ordinates. The non- vanishing components of the metric 
tensor are g rr = 1, g^ = r 2 sin 2 (#), and gee = r 2 . The 
Christoffcl symbols can be calculated from the metric 
tensor by using standard differential geometry relations. 
Note that our simulation region does not include the 
poles because there, the determinant of the metric tensor 
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FIG. 6. Torque coefficient as a function of the Reynolds num- 
ber for the case of two concentric rotating spheres. The the- 
oretical result has been taken from Ref. [lq] , a nd the experi- 
mental data have been collected from Ref. Ilq. [1711. 



becomes zero and therefore it is not possible to calculate 
its inverse. To circumvent this problem, we simulate the 
region 9 6 (tt/6, 5tt/6). We set r = 0.8 and use a lat- 
tice of size 32 x 1 x 384. In order to vary the Reynolds 
number, we change the azimuthal velocity d<p/dt. The 
boundary conditions have been chosen periodic for the 
case of 4>, and fixed for the case of r and 9. In the inset 
(left) of Fig. [2] we show the critical Reynolds number for 
different configurations which is in good agreement with 
the experimental values given in Ref. 0. In this figure, 
we can also observe the radial and polar components of 
the velocity, and see that there are two small vortices 
located at the equator and two large ones at high and 
low latitudes, in agreement with experiments and other 
numerical simulations [(| 0, E3] ■ 

We have also measured the torque coefficient defined 
by, 



T r = 27ra J 



er r sin 2 (9)d9 , 



(28) 



where a rc j, is the shear stress tensor, which in the context 
of lattice kinetic theory can be calculated by, 



(29) 



The torque coefficient is then computed via the follow- 
ing relation [l7j 



> 5 (f)' 



(30) 



find good agreement with the experiments. The small 
discrepancy can be due to the approximation taken in 
Eq. (|2"9"|) and the implementation of the boundary condi- 
tion. 

In order to study the Taylor- Couette instability for the 
case of two concentric rotating tori, which to our knowl- 
edge has never been done before, we use a lattice of size 
64 x f28 x 64 cells in the orthogonal coordinate system 
of the torus, (r, u, v), being r the radial, u the axial, and 
v the tangential coordinates. The Christoffel symbols 
and the components of the metric tensor can be read- 
ily calculated from differential geometry relations. The 
major radius of the tori has been taken as 4.0, in nu- 
merical units. The other parameters are the same as 
in the previous simulations, and to vary the Reynolds 
number we change the tangential velocity dv/dt. In this 
case, a and b are the minor radii of the inner and outer 
tori, respectively. We use periodic boundary conditions 
for the coordinates u and v, and fixed boundaries for r. 
In addition, the critical Reynolds numbers for different 
configurations can be observed in Fig. [21 showing values 
around 10% larger than for the case of cylinders. 



In Fig. |nj we show the comparison between our re- 
sults, the theory for Re — > 0, and the experiments. We 
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